https://bookdown.org/mcwimberly/gdswr-book/raster-data-discrete-variables.html

In this R practice, we will work with importing and managing the National Land Cover Dataset (2015).

Packages needed (download if you are missing any)

Next, let’s import some of the files we will need to work with.

NLCD <- raster("~/Desktop/nlcd/kad_NLCD_4-7/nlcd_kad_37.tif")
kad_data <- read.csv("~/Desktop/landscape_genetics_2022/labs/landgen_spatiallab_week3/kad_localities.csv")
kad_watershed<- sf::read_sf(("~/Desktop/aculley_KAD_master_folder/KAD_working_shapefiles/watersheds.shp"))

NLCD_shape <- sf::read_sf(("~/Desktop/nlcd/kad_nlcd_test4-7/nlcd_kad_test.shp"))
plot(NLCD)

rasterdf <-function(x, aggregate = 1) {
 resampleFactor <- aggregate        
inputRaster <- x    
 inCols <- ncol(inputRaster)
inRows <- nrow(inputRaster)
 # Compute numbers of columns and rows in the new raster for mapping
 resampledRaster <- raster(ncol=(inCols / resampleFactor), 
                       nrow=(inRows / resampleFactor))
  # Match to the extent of the original raster
 extent(resampledRaster) <- extent(inputRaster)
 # Resample data on the new raster
 y <- resample(inputRaster,resampledRaster,method='ngb')

  # Extract cell coordinates  
 coords <- xyFromCell(y, seq_len(ncell(y)))
 dat <- stack(as.data.frame(getValues(y)))
 # Add names - 'value' for data, 'variable' to indicate different raster layers
 # in a stack
 names(dat) <- c('value', 'variable')
dat <- cbind(coords, dat)
 dat
 }
nlcd16_df <- rasterdf(NLCD)

ggplot(data = nlcd16_df) +
  geom_raster(aes(x = x, y = y, fill = value))

This map and the accompanying legend are not particularly useful. The values stored in the land cover raster are numeric codes. For example, a value of 11 represents water wheres a value of 21 represents low-density residential development. The default chart created by ggplot() maps these values on a continuous scale rather than as categories, and the resulting symbology and legend have no real meaning. A detailed list of all the NLCD codes can be found at the following link: https://www.mrlc.gov/data/legends/national-land-cover-database-2011-nlcd2011-legend.

To generate a better-looking map, will start by creating two vectors - one containing the numeric codes for all 16 NLCD classes that occur in the conterminous United States, and another containing a human-readable names in the same order as the codes.

LCcodes <- seq(1,19,1)

LCnames <-c(
  "TemperateSub-polarNeedleleafForest",
  "Sub-polarTaigaNeedleleafForest",
  "TropicalSub-tropicalBroadleafEvergreenForest",
  "TropicalSub-tropicalBroadleafDeciduousForest",
  "TemperateSub-polarBroadleafDeciduousForest",
  "MixedForest",
  "TropicalSub-tropicalShrubland",
  "TemperateSub-polarShrubland",
  "TropicalSub-tropicalGrassland",
  "TemperateSub-polarGrassland",
  "Sub-polarPolarShrubland-lichen-moss",
  "Sub-polarPolarGrassland-lichen-moss",
  "Sub-polarPolarBarren-lichen-moss",
  "Wetland",
  "Cropland",
  "BarrenLands",
  "Urban",
  "Water",
  "SnowAndIce")


LCcolors <-c("#009300",
  "#949C70",
  "#006300",
  "#006300",
  "#148C3D",
  "#5C752B",
  "#B39E2B",
  "#B38A33",
  "#E8DC5E",
  "#E1CF8A",
  "#9C7554",
  "#BAD48F",
  "#408A70",
  "#6BA38A",
  "#E6AE66",
  "#A8ABAE",
  "#DC2126",
  "#4C70A3",
  "#E1FAFF")
names(LCcolors) <- as.character(LCcodes)
LCcolors
##         1         2         3         4         5         6         7         8 
## "#009300" "#949C70" "#006300" "#006300" "#148C3D" "#5C752B" "#B39E2B" "#B38A33" 
##         9        10        11        12        13        14        15        16 
## "#E8DC5E" "#E1CF8A" "#9C7554" "#BAD48F" "#408A70" "#6BA38A" "#E6AE66" "#A8ABAE" 
##        17        18        19 
## "#DC2126" "#4C70A3" "#E1FAFF"
#LCcolors <- attr(NLCD, "legend")@colortable[LCcodes + 1]
#names(LCcolors) <- as.character(LCcodes)
#LCcolors
NLCD_df <- fortify(NLCD_shape)

mapping shapefile

ggplot(data = NLCD_shape) +
  geom_sf(aes(fill = gridcode), color = NA)

Now the ggplot() function can be used to plot the NLCD land cover raster. In the aes() function, the raster value is converted to a character so that ggplot() will recognize is as a categorical variable. Thescale_fill_manual() function is then used to specify the colors that match the NLCD codes. We are removing the second item of LCnames because there are no “Perennial Ice/Snow” areas in Georgia, so we do not want this class to show up in the legend. (You can check this out by running unique(nlcd11) and noting that “12” is not present.) The theme() function is used to remove the axis titles and change the appearance of the background.

ggplot(data = nlcd16_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = LCcolors,
                    labels = LCnames,
                    na.translate = FALSE) +
  coord_sf(expand = F)